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Abstract 



The aim of this paper is to develop a mathematical framework for opto-elastography. 
In opto-elastography, a mechanical perturbation of the medium produces a decorrelation 
of optical speckle patterns due to the displacements of optical scatterers. To model this, 
we consider two optically random media, with the second medium obtained by shifting the 
first medium in some local region. We derive the radiative transfer equation for the cross- 
correlation of the wave fields in the media. Then we derive its diffusion approximation. 
In both the radiative transfer and the diffusion regimes, we relate the correlation of 
speckle patterns to the solutions of the radiative transfer and the diffusion equations. We 
present numerical simulations based on our model which are in agreement with recent 
experimental measurements. 
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q ; 1 Introduction 

When a strongly scattering medium is illuminated by a coherent laser, the transmitted light 
forms a speckle pattern whose statistical properties depends on the statistical properties of 
| the medium. When the optical scatterers are displaced in some region, the speckle pattern 

changes and its correlation with the original speckle pattern depends on the amplitude of the 
displacement and on the local optical properties (scattering and absorption) of the region 
affected by the displacement. 

Focused ultrasound introduces such displacement of scatterers through two different mech- 
anisms. First, ultrasound focusing generates oscillating compressive strain in the focal region; 
the oscillation of the optical scatterers is in the MHz frequency range. Second, high-intensity 
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focused ultrasound may also generate low frequency (of kHz range) strain in elastic media, 
which in turn generates shear wave propagating in the medium. Both modifications of scatter- 
ers can cause decorrelation of speckle patterns. However, when speckle patterns are recorded 
with an exposure time sufficiently long with respect to the compressive oscillation while re- 
mains short with respect to the shear motion, the effect of the high frequency compressive 
motion averages out and the decorrelation in speckle patterns is dominated by the second 
mechanism. Based on this idea, transient opto-elastography experiments have been carried 
out by Bossy et al. in [HE], where they found qualitative relations between the decorrelation 
of speckle patterns and optical absorption as well as mechanical properties (such as Young's 
modulus) of soft biological tissues. 

The main objective of this paper is to provide an analytical model that relates the decor- 
relation of speckle patterns to the displacement of the scatterers. Such a model hence helps us 
to understand the aforementioned " acousto-elasto-optic" phenomenon. To this end, we start 
with the Helmholtz equation with the refractive index having a highly oscillatory random 
part. We consider the regime where the correlation length of the random medium is of the 
same order as the wavelength and both are smaller than the typical propagation distance. 
We denote by e the ratio between the correlation length of the medium and the propagation 
distance and assume therefore that £<1. We also assume that the relative amplitude of the 
random fluctuations is weak, of order yfe. This is known to be a scaling regime where the 
random medium interacts with the propagating high-frequency waves [11] . We consider two 
random media, with the second medium obtained by shifting the first medium in some local 
region. 

The correlation of wave fields is well described by the Wigner distributions [9] . Following 
the techniques of [IT], we formally derive the radiative transfer equations (RTEs) for such 
cross-correlations of wave fields in two random media in the limit as e goes to zero. Radiative 
transfer limits for waves in two different random media have been considered, e.g. in [U 12], 
but the case of two media related by a local shift considered in this paper is new. The 
salient effects of this local shift of random media are: it introduces a phase modulation to the 
scattering kernel of the RTE in the region affected by the shift; further, when the amplitude 
of the shift is large, say much larger than the wavelength, by non-stationary phase the RTE 
for the cross-correlation function is intrinsically absorbing in this region. Next, following the 
techniques of [8] O [7] , we derive the diffusion approximation of the RTE in the regime when 
the mean free path is small; this simplification is useful for numerical simulations. In the 
special case of large shift, the cross correlation vanishes in the region affected by the shift. 
We shall see that this accounts for the loss of correlation of the speckle patterns. 

The paper is organized as follows. In the next section, we present the model of two 
random media considered in this paper; they are related by a local shift. In Section [3l we 
derive the RTEs for the Wigner distributions of the wave fields in the aforementioned two 
random media, in the limit when e goes to zero; we also present the formula for speckle 
pattern correlation in terms of the solutions of these RTEs. In Section H] we derive the 
diffusion approximation of the RTEs in the limit when the mean free path goes to zero, and 
again we derive the formula for the speckle pattern correlation in the diffusion regime. In the 
special case of large shift, using the derived diffusion approximation and correlation formula, 
we run numerical simulations which confirm that the diffusion equation model is able to 
capture the loss of correlation of the speckle patterns. 
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2 Wave equation and heterogeneous media 



In the microscopic scale, light propagation is described by the Maxwell equations. In the case 
when scalar approximation is valid, these reduce to the following Helmholtz equation for the 
electric field u, 

A«(x) + /c 2 n 2 (x)u(x) = 0, (1) 

where ko is the wave number of the light in vacuum and the spatially varying refractive index 
n 2 (x) models the heterogeneous medium. 

The following model for n 2 is adopted: 

n 2 (x):=n 2 (l + 2aV(j)), (2) 

where V(x) is a real valued mean-zero random process. Therefore, the heterogeneous medium 
has mean refractive index n 2 ,, which is assumed to be a constant, and its relative fluctuations 
are captured by 2<rV(x/Z) where the numbers a and I model the strength and the correlation 
length of the fluctuations respectively. 

The random process V(x) is assumed to be stationary, i.e., statistically homogeneous. 
The two-point correlation function of this process is 

fl(x) = E[F(y)V(y + x)] = E[V(0)V(x)]. (3) 

where E stands for the expectation with respect to the distribution of the random medium. 
Throughout this paper, we adopt the convention that the Fourier transform of a function / 
on M. d is defined by 

Using this notation and the stationarity of V, one easily verifies that 

E[F(p)y(q)] = £(p)5(p + q), (5) 

where S is the Dirac distribution. From now on the product k 2 riQ is denoted by k 2 . 

Let L denote the typical propagation distance. In the high-frequency regime, the ratio 
between the wavelength and the propagation distance, denoted by e := k^ 1 /L, is much 
smaller than one. Let k := l/k -1 be the ratio of the correlation length of the random medium 
to the wavelength. We consider the so-called weak coupling high-frequency medium which 
corresponds to k = 1 and a = yfe. This is known to be a situation where the heterogeneous 
medium interacts with the high-frequency waves . Denoting by u £ (x) the scaled function 
u(Lx), the Helmholtz equation becomes 

yA W £ (x) + ^n £ (x) + v^y £ (x)u £ (x) = 0. (6) 

with V £ (x) = V(x/e). 

Two random media. As mentioned in the Introduction, we are interested in the cor- 
relation of two wave fields in two random media. With the application to opto-elastography 
in mind, we view these two media as configurations of scatterers at two instants of the shear 
motion. Let u\ be the first wave field which solves with Vf(x.) = V(x/e). Let u\ be the 
second wave field which solves ([6]) with V^ e (x) given by 

Vi (x) = Vf(x + £ 0(x)) = V{- + 0(x)), vm = V{-). (7) 
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Here <f> is a continuous and compactly supported vector field. V 2 e * s obtained from V{ by the 
application of the diffeomorphism x — > x+e0(x). Hence 4> can be thought as the displacement 
field that shifts the configurations of optical scatterers. Note that the amplitude of the shift 
is of order e, which is of the order of the optical wavelength. 

Correlation of waves. It is well known that the correlation of two scalar wave fields u 
and v is well described by the Wigner transform IU[-u,-u], which is defined by 

W[ti,«](x,k) :=(^p/ e 4k 'Mx-f )v(x+^)dy, (8) 

where the bar denotes the complex conjugate. According to (j3J), the Wigner transform may 
be thought as the Fourier transform of the two-point correlation function of the fields. In 
this paper we are interested in 

W£(x,k) = W^ } «f](x,k), (9) 

where the wave fields id, j = 1,2, are described earlier and they correspond to two random 
media Vj related by ([7]). We note that 

J W7,(x,k)dk = uf(x)5f(x). (10) 

In particular, fW^dk is the energy density |uf| 2 and W{ x itself corresponds to frequency- 
specified energy density of the first wave field. The functions W 22 and Wf 2 can be interpreted 
similarly. 

In the next section, we derive the radiative transfer equation (RTE) for in the limit 
e — > 0, and we relate the correlation of speckle patterns to these Wigner distributions. 



3 Radiative transport equation for the Wigner distributions 

The goal of this section is to derive the RTE for the correlation of waves. Recall that the two 
wave fields u|, j = 1,2, solve 



e 2 1 

-A^(x) + -^ £ (x) + v^y/(x)^(x) = 0, 



(11) 



where V 2 an d satisfy ([7]). Recall that is defined in ([9]). For the moment, we assume 
that u £ j has prescribed plane wave behavior at infinity to simplify the calculation and derive 
RTEs for the limits of WJ^'s. For bounded domains, these RTEs are valid in the interior of 
the domain, and appropriate boundary conditions should be imposed. 

Using ([9]) and (fTT|) . integration by parts and the prescribed plane wave behaviors of Uj's, 
we find 
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Formal derivation of RTE for W{± is classic; see We observe, however, that the 

equation for Wf 2 and W 22 are not standard because of the shift 4>. In fact, correlation of 
wave fields in different random media has been considered in [IJ, but they do not cover the 
case when one medium is a local shift of the other as described in (|7|) . 

Now we adapt the formal derivations in [T] to derive RTE for the cross-correlations 
W[ 2 and W 22 . In the process of the derivation and as in the references, various assumptions 
will be made whose rigorous justifications are known to be difficult. To form a closed equation 
for Wf 2 , we rewrite the equation it satisfies as follows using the Fourier transform of V, 



ik-VW£j(x,k) = -= 



dydp 



(2vr) c 



V(p) 



ip-d+l+^x+f )) _ e -ip-(f-i) 



uf(x-^)u|(x+^). 



We replace the function <^(x + ey/2) in the phase function by 0(x) and neglect the error 
term of order e. It follows that 



flc-VWf 2 (x,k)«-^ 
1 



dydp 



(2ir) d 

-p-ty(p) 



v(p) 



ip.(| + i+0(x)) _ e -tp-(f-| ) 



nf(x-^)n|(x + ^) 



e -^«A(x)^ 2 ( X)k _ £) _ ^ £ 2 (x,k + P 



dp, 



where in the second equality we have used the definition 
have 



2' 1ZV 2 
JJ). By the same argument we also 



ik-VW| 2 (x,k) = 4= / e- ip -(? + ^ (x) )y(p 
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dp. 



3.1 Radiative transfer limit by multiscale expansion 



In this subsection we consider the limit of W( 2 and W 22 as e goes to zero using a formal 
multiscale expansion argument. The equations derived earlier can be written as 



k • VW£(x, k) + -Lp^fo k) = 0. 



(12) 



Here jl takes the value 12 or 22. The operator Vji is defined accordingly by 

V 12 Wf 2 = i J e-^V(p) [e-*P'*Ww^(x J k-|)-Wf 2 (x,k+|)]dp, 
VnW& = % / e- i P-(f + ^ x ))y(p) K £ 2 (x,k-?)-^| 2 (x,k + ?)l dp. 



Following the method of multiscale expansion, we define the fast variable £ = x/e, and 
assume that the following expansion is valid: 



Wl (x, k) = Wf (x, £ , k) + yfeWf (x, £, k) + eiyjf } (x, £ k) + 



(13) 



As e goes to zero, the leading-order term dominates, and the goal is to derive a RTE 
for this term. To this end, we substitute the multiscale ansatz (|13|) and the relation 

v = v x + -v € 

s 
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into (|12p . Equating terms that are of equal order in e, we find: 

o(e- 1 ) ■. k- V/:ir;;" - 0. 

0, 



k-V € < 



0(1) 



k • V x ^ 0) + k • VjWf + 



0. 



(14) 
(15) 
(16) 



The strategy of the derivation is as follows: The first equation imposes that is inde- 

pendent of the fast variable £. The second equation can be explicitly inverted and provides 
a representation of Wjf in terms of Wjf . Take statistical averaging on the third equation, 
and observe that 

E[k-V^ 2) ] =0, 

which is true because the statistics of V is homogeneous in the fast variable £; see the 
statement above ([3]). We also note that it is natural to assume that W^f is not random. 
Then the third equation reduces to 

k-V x W$ ) +E[P jl wU}=0. 



Therefore, it suffices to invert (|15p and evaluate the expectation above. 

Limiting equation for Wf 2 - Upon adding an absorption term 0W^ for regularization, 
we rewrite ([TBI as 



^(x) W rW (k _£ ) _ w g) Ck+ E ; 



dp = 0. 

Taking Fourier transform in the fast variable £, we obtain from the above equation that 

V(C) \e-«-*MW®(k - §) - W$\k + | 

r(i-) ■ - ' 
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k • C + it 



(17) 



To estimate the expectation of VizW^ > we rewrite this integral using the Fourier transform 
of and get 
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Using (|17p . we can represent this integral in terms of as follows: 
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dpdq. 
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Taking expectation, using the definition ([5]) and the fact that R(p) = R(—p), we obtain that 
the expectation V^W^ can be evaluated as 



(k) - e-^M w$ (k - p) e'P^M ) (k + p) - W$ (k) 



* / * (p) j ^ -(k- f) .p^ " — d » 

1 ^ -i(|k| 2 -|p| 2 ) + ^ lP ^ -i(|k| 2 -|p| 2 )-^ P 



0->O 



> 4vr J R{p - k)[W$\k) - ^ 1 ( 2 ) (p)e i(p " k) ^ (x) ]5(|k| 2 - \p\ 2 )dp. 



In the first equality, we have used the change of variables k — p i— > p and k + p i— > p for the 
two integrands respectively. In the last equality, we have used the fact that 

-> ird(x) 



x 2 + 9 2 



as a distribution of the one-dimensional variable x. Finally, we obtain the following RTE for 



k- Vw£>(x,k) +4ir J R{p-k)[W^\k) - iy i ( 2 ) (p)e i ( p - k )^W]5(|k| 2 - |p| 2 )dp = 0. (18) 

Limiting equation for W^- The above procedure can be applied to the equation 
satisfied by W^i as well. Upon adding a regularizing term and using Fourier transform, we 
solve equation (fT5j) by 

W-gW.k) = 1 r t , (19) 



k • C + it 

Using this solution and the Fourier transform representation, we find that 7^22^22^ nas ^ ne 
form 



7 ? 22wg ) (x,€ > k) = i ^ e- i ( p+q )-^+^ x »y( P )y(q) 



w 2 ( 2 0) (k-f-§)-w 2 ( 2 0) (k-f + §; 



(k-§)-q+i0 

iy 2 ( 2 0) (k+£-a)-^ 2 °)(k+£ + a) 

(k+f)-q+ifl 



(ipciq. 



This expression is much simpler compared with that of V^Wy^ because the phase modifica- 
tion due to <p is uniform for the W^) above. Taking expectation and using ([5]) we find that 
this modification has no effect on the expectation. In fact, we have 

E[P 22 ^ 2 ( 2 1} (x,^,k)] -> 4vr y R(p - k)[I^ 2 ( 2 0) (k) - w£\p)]5(\k\ 2 - |p| 2 )dp. 
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Finally, we obtain the RTE for W$: 

k ■ W<°>(*, k) + 4. / *<p - k )K»(k) - <<p)WW - IPI'WP = 0. (20) 

We note that this is the same as the classic RTE limit for wffl . 

Summary of the results. Now we summarize the above results and discuss the bound- 
ary conditions when the problem is posed on a bounded domain as often encountered in 
practice. Noticing that 

« 2 -| P | 2 ) = ^(|k|-|p|), 
we define the scattering kernel 

(j(p,k;|k|) = 2vr|k| d - 3 i?(p-k), (21) 

and the absorption coefficient 

S(k) = ^ J 5(\k\ 2 - \p\ 2 )R( P - k)dp = J «r(p, k; |k|)dp. (22) 

Here and in the sequel, the vector p is defined to be the direction p/|p| for non-zero p. 
Abusing notations we use Wj\ for the leading-order terms Wjf . With these notations, the 
equations (fT8|) and (p0|) become 

k-VW 12 + ZW 12 = [ ( r(p,k;|k|)e i l k l(P- lc ^( x % 1 2(x,|k|p)dp, (23) 

k • VWjj + EWtf = / <t(p, k; |k|)^-(x, |k|p)#, j = 1, 2. (24) 



On a bounded domain X C M. d , we need to equip the RTEs with proper boundary 
conditions. We consider the monokinetic case when |k| is fixed, so the energy density is 
specified by the spatial variable x and the direction k. The RTEs are therefore posed on the 
phase space X x S d ~ l , with |k| as a fixed parameter. We define the incoming boundary T_ 
and the outgoing boundary T + as 

r± := {(x, k) | x € dX, ±k • i/(x) > 0}, (25) 

where is the outer-pointing normal at the point x on the boundary dX. When a 

bounded domain is considered, the RTEs ([23]) and (|2"3]l for the Wigner transforms Wji's 
should be understood as for (x, k) G X x S^ 1 with the boundary condition 

W i ,(x,k)=p(x,k) ) (x,k)eL, (26) 

where p models the light intensity incoming at point x € dX in the direction k. In the case 
of an incident laser beam the support of p is spatially limited to x € dXi where dXi is the 
part of dX where the laser beam is applied. 

The case of large shifts. When the shift 4> is much larger than the wavelength, i.e., 

|k||0|»l, (27) 
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the RTE for W± 2 can be simplified. Indeed, since the integral is taken over the sphere which 
has non- vanishing Gaussian curvature, we can apply [12\ Theorem 1.2.1] which may be viewed 
as an analog of the Riemann-Lebesgue lemma, and conclude that the integral on the right- 
hand side of (fTSj) is of order (|k| |0|)~( d_1 )/ 2 and hence approaches zero. Consequently, (|18|) 
should be modified as follows 

k • VW12 + T,Wi2 = 0, (x,k)a s xS" 

f j 9 (28) 

k-VW 12 + ZW 12 = / <7(p,k;|k|)Wi 2 (x,|k|p)dp, (x, k) 6 X£ x S . 

Here X s is the support of 0, i.e., the region in which the scatterers are shifted, and X£ is the 
complement of X s in X. 



3.2 Speckle pattern correlations in the RTE regime 

As we have seen, the Wigner transform W[«,if] of a wave field u has the interpretation of 
direction-resolved energy density. Hence RTE for W[u, u] is a very good model for light 
propagation. In this section, we relate the correlation of optical speckle patterns to integrals 
of Wigner transforms. 

In the opto-elastography experiment, emitted light intensity is measured at a part dX m 
of the domain boundary dX and the data are {|u^(x)| 2 | x G dX m }, j = 1, 2. The correlation 
of two speckle patterns u\ and u\ is defined by 



C 



12 



m\ 2 - {\av)){wiv - {\uim 



vm\ 2 -(K\ 2 )) 2 )vmF 



(Kl 2 )) 2 ) 

where {A) denotes the spatial average over the boundary dX m , that is 



(A) 



\dX n 



A(x)dx. 



(29) 



(30) 



Note that in the above equation, the notation dx means the induced Lebesgue measure on 
the boundary dX m , and |9X m | is the area of the boundary. 

Assume that the complex amplitudes (u\ C 2 -valued random process, satisfy the 

circular symmetric Gaussian distribution, and assume also that the spatial average can be 
thought as ensemble averages (taking expectations), we have 

ei2\ /l„.ei2\ 



(\u^\u^) = (uy i )(u^) + (\u c /)(\u c l 
Using this equality, we find that the numerator of (j29[) is 

2 



3,1 = 1,2- 



(31) 



K«!«§>l' 



1 



\0X n 



u\ (x)«2(x)dx 



1 



\dX„ 



W 12 (x,~k)dkdx 



where r mj+ := {(x, k) | x G 3X m ,k • u{x) > 0}. In the second equality above, we have used 
the fact that u\u\ is the integral over k of W\ 2 {x, k), as seen in (fT0|) . Since the product u\ «| 
only accounts for the outgoing light, we take the integral of W\ 2 over outgoing directions. 
Similarly, the denominator in (j29j) is given by 



<i«!rxi«s 



ei2\ 



1 



W u (x,k)dkdx 



W 22 (x,k)dkdx. 
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Combining these calculations above, we find that in the RTE regime, the correlation of two 
speckle patterns is: 





J 


( W i2 (x 


k)dkcbc 


2 


f Wu(x,i)dkdx I 





We remark that modeling of the outgoing light depends on the measurement set-up. In 
the above model, outgoing light in any directions are completely captured. If light is collected 
in the far field by a lens, then we should only take the integral in k in a cone that corresponds 
to the aperture of the lens. Such a situation does not affect qualitatively the results in the 
sequel. 



4 The diffusion limit for the radiative transfer equations 

The RTEs derived in the previous section are posed on the phase space X x S rrf-1 . This 
causes difficulties for instance for numerical simulations. In this section, we consider their 
diffusion approximations which are much easier to deal with for the purpose of simulations. 

For the sake of simplicity, we consider the case when the scattering kernel a, as defined 
in (|2ip . takes the form 

o"(p,k; |k|) = cr(p • k; |k|). 

In this case, the total absorption section £ defined in (|22|) depends only on the parameter 
|k|. In the sequel of this paper, we will use the normalized scattering kernel 

/(p ■ k; |k|) = — j— a(p-k;|k|). (33) 

Very often and as we have already done so, the dependence on |k| is not written explicitly. 
The diffusion approximation is valid in the regime where the mean free path rj = 1/S is much 
smaller than the size of the domain. Here we assume r\ <C 1 and the domain is of order on^\. 

In the following two subsections, we adapt the classical derivation of diffusion limits of 
RTE [8, 3, 7J to the case considered in this paper. The point here is that we have to manage 
the shift field cf>. Though rigorous derivation of diffusion limit is possible, we do not pursue it 
here and adopt the formal multiscale expansion argument only. We consider two interesting 
situations. In the first one, the amplitude of the shift cf> is very small so that |k||0| is of order 
ij. In the second one, the amplitude of the shift cf> is very large so that |k||0| is much larger 
than one. As seen in the last section, this leads to the RTE ()28|) . 

4.1 The case of small shift 

In this section, the amplitude of the shift 4> is assumed to be small so that |k|0(x) = r]tp(x.). 
The diffusion limit for W\\ and W22 are classic [3 Section XXI. 5. 4]; hence we concentrate on 

Alternatively, when rj is bounded but not necessarily small, diffusion limit is still valid when we consider 
the problem on a very large scale and rescale the spatial variable. 
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that for W12 in (|23j) . Dividing on both sides of this equation by S, we get 



r/k- VWi 2 (x,k) 



S d-i 



/(p • k) [ e ii(p-k)*)^ 12 ( X) ^) _ Wi 2 (x,k) 



dp. 



We follow the idea used in [T] and define a new function which takes account the phase shift 

Wi 2 (x,k) = e^ (x) Tyi 2 (x,k). 
Then it is easy to verify that W\2 satisfies 



r/k • VW 12 - «r/ 2 S a ^i2 



S d-i 



/(p-k) W 12 (x, P )-W 12 (x,k) 



dp. 



(34) 



Here, the intrinsic "complex absorption" coefficient is given by T, a = Tr(k ® kV?/>(x)), where 
k k is the projection matrix kk* and Tr means taking the trace. 

Since W\2 approximates W\2 as 77 goes to zero, it sumcies to consider the limit of TUi2- 
Abusing notations, we still denote this function by W\2- To start the formal derivation, we 
substitute the ansatz 

into the equation (|34|) . Equating the terms that are of equal order in 77 we find 



(35) 



0(1) 







/(p-k) w£\ P )-w£>(k) dp, 



(0), 



0(r/): k-Vw} 0) 



S d-i 

.(o) 
12 



/(p-k) ^(^-^(k) 



dp, 



<D(rj 2 ): (k-VWjJ-ffi^ 



,(0) 
12 



/(p-k) w£\p)-W£(k) 



S d-i 



(2), 



dp. 



To study these equations, we define the integral operator 



Kh{k) 



S d-i 



/(p • k)/i(p)dp. 



(36) 
(37) 
(38) 

(39) 



Using this definition, Eq. (|3"B"j) can be recast as (K-I)W^ = 0. We verify that J sd ^ /dp = 1; 
we assume also that / is uniformly bounded from up and below by positive numbers. In this 
case, the constant function is known to be the only eigenvector of K corresponding to the 
eigenvalue one [7J. Furthermore, the integral equation (K — I)h = v is solvable only if 
f S d-i vdp = (Fredholm alternative). 

Using these results, Eq. (f36|) shows that does not depend on the direction variable 

k, so W$ = W^\x.). Eq. (f37|) relates wj^ to W 12 . In fact, we need to solve the integral 
equation 



(K-I)W$=k.VW^( X ). 
Let hj(k) be the unique solution (whose integral over S^ -1 is zero) to the equation 

(K - I)hj(k) = ej -k, (40) 

where ej is the j'-th unit vector in the orthonormal basis of W 1 . Indeed, this is possible 
because ej ■ k integrates to zero on the sphere S^ 1 . Let h(k) be the vector field (hi, /i 2 , ^3)*- 



It follows that W[l\x,k) = h(k) • VlU^(x). 



(o). 
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Substituting this representation of into Eq. (|38|) and integrating the equation over 
,S d_1 , we find that the right-hand side vanishes due to the symmetry of / and we get 

/ k- V [h(k) • Vl¥ (0) (x)| dk-i [ Tr(k^kV'0(x))iy (o) (x)dk = O. 
J S d-i L J Jgd-i 

Carrying out these integrals, we find 

- D : VV^ (0) (x) - ^(V • </>(x))ty (0) (x) = 0. (41) 

Here, the symbol : denotes the Frobenius inner product of two matrices, d is the space 
dimension, and w,i is the area of the sphere S d ~ x . The diffusion matrix D is given by 



D := — h(k) <g> kcik. (42) 
To complete the derivation, it suffices to find hj(k) that solves (|3D)>. Define 

^-T^ilfey (43) 

with ^ 

</(|k|) = / /(p • e; |k|)p • edp = C d [ /(//; |k|)Mu, (44) 
Js*- 1 J-i 

where e is a unit vector in S^ -1 and Cd is a constant depending only on the dimension. We 
note that g is independent of the choice of e and that g < 1. Let Qj be an orthonormal 
matrix so that Qjk = Gj. We check that 

(K - J)(-k • e,-) = k • a, - / /(p • k)p • e,-dp 

= k • ej - / f(Qjt> ■ ej)Qjf> ■ QjBjdp = k • — / /(p • ey)p • Qjejdp. 

In the last equality, we changed the variable Qjp to p. We have the decomposition 

QjQj — [gj * QjGj^Gj -f - CGjj_ — [Qjk • QjGj^Gj -\- CGjj_ — (k • GjjGj -\- CGjj_ : 

where cGj± is perpendicular to Gj. By symmetry, the contribution of cej± to the spherical 
integral above vanishes. We can then check that 

(K - I)(-k • Gj) = k • gj - (k ■ ej) / f(p ■ Gj)i> ■ ejdp = (1 - ^(|k|))ej • k. 
Since /t,- defined in (|43p clearly integrates to zero on 5' d_1 , the above calculation shows that 



hj solves (}4T)|) . Using this solution for D in (|4"2j) . we find that D = zudl/[d(l — g(|k|))] where 
I is the identity matrix. In dimension three, the limiting diffusion equation for which 
is denoted by Wu becomes 

- V • 3(1 _ 4 J (|k|)) W 12 (x, |k|) - *y (V • V(x))Wi2 = 0. (45) 
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For the autocorrelation functions W\\ and W22, we take ij) above to be zero and recover 
the classic limiting diffusion equation 



v - 3 (i- 4 9 Vi)) w " (x - |k " =o - j=i ' 2 - < 46 » 



The above diffusion equations (|45j) and (j46]) should be equipped with proper boundary con- 
ditions. We impose that 

Wji(x) = g(x), x G dX. (47) 

Here, g(x) models the incoming light intensity at the boundary. It can be derived from the 
boundary condition p(x, k) in (|26|) of the RTE. If p is independent of k, then q = p. The case 
when p is anisotropic requires a careful boundary layer analysis as developed in [HE], which 
shows that there exists a linear operator (local in x) that maps p to q. In the case when / is 
constant (isotropic scattering) and p(x,k) = p(x, //) where fi = — k • i'(x), this map is given 
by 

?(x)= / p(x, M )ff(M)^A*, (48) 
J 2 

where is the so-called Chandrasekhar //-function. The evaluation of H can be found, 

for instance, in Section 1.5]. As a result the support of q is spatially limited to the part 
dXi of the boundary dX where the laser beam is applied. 



4.2 The case of large shifts 

In this section, we consider the more practical case when the shift cj) of the scatterers has 
large amplitude so that |k||^»| 3> 1. As we have derived in Section 1331 the RTE for W\ 2 takes 
the form (|28p . Consider the limit 77 = 1/S <C 1, these equations can be written as 



r]k ■ VW12 + W12 = 0, (x, k)£l s x S' 



d-X 



r/k • VW 12 + W 12 = [ /(k, P )^ 12 (x, |k|p)dp, (x, k) e A s c x 



On the unshifted region Xg, the second line of the RTEs above takes the classic form 
and its diffusion limit is (|46p in three dimensions. On the other hand, sending r] to zero, we 
deduce from the first line of the equations above that W12 goes to zero in the shifted region 
X s in the limit. Consequently, the diffusion limit for the cross-correlation W\ 2 in the case of 
large shift is: 

- v 'I(T^W™ k l> = °- (49) 

W 12 = 0, x G A s , 

and it is equipped with the boundary condition (|47p . The diffusion limits for W\\ and W22 
remain unchanged; that is, they are given by (|46p with boundary condition (|47p. 



4.3 Speckle pattern correlation in the diffusion regime 

In this subsection, we revisit the formula (|32p and rewrite it in terms of the functions in the 
diffusion limit. As before, we denote by dX m the part of the boundary where measurement 
is taken. We naturally assume that the part dXi of the boundary that is illuminated by the 
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incident laser beam has an empty intersection with the part dX m where the outgoing light 
is measured. Therefore Eq. (|47p implies that the leading-order term Wji(x) in the diffusion 
approximation of the direction-resolved correlation function Wji(x, k) vanishes on dX m . It 
is then necessary to look for the first-order corrections in the expansion of Wji(x, k). This is 
discussed in detail in [10] and we follow the ideas there. 

In the diffusion regime, using the expansion (|35p with found in section 14. 1[ we can 
write Wji (x, k) as 

Wjfa k) = w^,(x) - 1 _ 71 gm * ■ VW^-,(x) + • • • , (50) 

where Wji (x) is the function involved in the diffusion equations (|4"5j) (06]) fi9]) . This expansion 
is valid only inside the domain X and not at the boundary dX, again because of the presence 
of a boundary layer which gives rise to a correction of order ij that cancels the first-order 
corrective term in (|50p for x 6 dX m and i/(x)-k < 0. If, however, we carry out the calculation 
with this expansion formally, then we find that 



L 



I^(x,k)dk = - I f kdk • V^(x) = -— ^j— i/(x) • V%(x) 



k-i/(x)>0 1 - A-i/(x)>0 

where is a constant that depends on the dimension. The second equality follows from the 
decomposition k(k • i/(x))i/(x) + kj_ where kj_ is perpendicular to v, and the fact that the 
contribution of kj_ averages to zero because of symmetry. 

In fact the result obtained in this formal way is essentially correct up to the value of 
the constant Cd that may depend on the form of / and can be found numerically when / is 
constant in particular |10[ Sec. X.A]. It follows that the correlation of speckle patterns in the 
diffusion regime is given by 



C 
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2 




[ u{x) ■ VWi 2 (x)dx 






JdX m 






VWn(x)dx f v(x) 


VW 22 (x)dx 


JdX m 


Jax m 





(51) 



We remark that in the above formula, the constant C^f]) (1— g) is cancelled out, so the desired 
result does not depend on its value. 



4.4 Numerical simulations of the diffusion equation model 

In this subsection, we show some numerical simulations which confirm that the diffusion 
equation model derived above for the cross correlation W\% captures the loss of correlation 
in speckle patterns. The numerical simulations are in accordance with the experimental 
measurements published in [U [5] . 

Let the domain X be a two-dimensional square (—1, 1) x (—1, 1), and let a sequence of 
circles S{r n ) centered at (0,0) with increasing radius r n model the wave front of the elastic 

(n) 

wave introduced by the ultrasound modulation. Let C]y be the correlation calculated as 
in (|5ip on the right side of the square domain X where the two random media are those 
when the wave fronts are at S(r n ) and S(r n -i) respectively. Since <fi models the shift of the 
scatterers of these two random media, the support of <£> is the union of the supports of the 
elastic waves at the two instants, and it is enclosed inside the circle S(r n ). 
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To evaluate C\ 2 > we need to calculate Wn,W22 and W\2- For the first two functions, 
we solve the diffusion equation (|46p with unit diffusion coefficient on the whole domain X 
with boundary condition (|47p . which is taken as q = 1 on the left side and q = on the 
three other sides of the square domain X. For W12, since it vanishes on the support of 4> 
which has outer boundary S(r n ), we solve the first equation of ()49j) on the exterior of the ball 
enclosed by S(r n ); the boundary condition is (|4"T|) on dX and W12 = on the inner boundary 
S(r n ). These configurations of computational domains and elastic wave fronts are illustrated 
in Fig. HJa). 




(a) (b) (c) 



Figure 1: Computational domains and elastic wave fronts, (a) The dash lines are the wave 
fronts 5'(r ra )'s. (b) An optical absorber with radius r a = 0.2 centered at (0, 0). (c) An optical 
absorber with radius r a located at (0,0.1). 

To demonstrate the effect of optical absorbers, we also consider the case when such an 
absorber with radius r a = 0.2 is located at (0,0) and (0,0.1) respectively, as illustrated in 
Fig- H]-(b)(c). They are referred to as the case with centered and non-centered absorbers 
respectively. In these cases, the equation for W\\ and W22 are solved outside the absorber 
because light is completely absorbed inside. For W12, the equation is solved outside the union 
of the absorber and the ball enclosed the circle S(r n ). 

Note that for simplicity, the diffusion coefficient D and the boundary condition q on the 
left boundary are chosen to have unit value. This does not affect the results of the simulation. 
Indeed, the diffusion constant is cancelled out in (|5ip ; further, when the boundary condition 
q in (|48|) is uniform in x on the left side, its constant value will be cancelled out in (|51|) as 
well. 

Finally we plot as a function of r n in Fig. [2j We see that when there is no optical 
absorber, the correlation of speckle patterns drops immediately when elastic wave forms. 
On the other hand, when optical absorber is present, the correlation starts to decay only 
after the elastic wave exits the absorber, that is at r = 0.2 for Fig. QJb) and at r = 0.1 for 
Fig. [He). Hence the decay of correlation is sensitive to the location of the absorber. This 
can be exploited further for medical imaging purposes. 
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